import sys
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
import timeit
import cv2 as cv
from skimage import exposure
from scipy.fftpack import fft, fftshift, ifft, fft2, ifft2 # can replace scipy.fft to older version called scipy.fftpack
import scipy
from scipy import signal
from scipy.signal import cont2discrete, lti, dlti, dstep
from numpy import cos,pi,real
from matplotlib.pyplot import subplot, figure, title, stem, axes, plot, imshow
from scipy.signal import butter , lfilter, freqz
from scipy.interpolate import CubicSpline
def MSE_SNR_PSNR(SignalA, SignalB):
mse = (np.square(SignalA - SignalB)).mean()
error_ = (SignalA - SignalB)
sig_mean = (np.square(SignalA)).mean()
noise_mean = (np.square(error_)).mean()
snr = 10*np.log10(sig_mean/noise_mean)
psnr = 20*np.log10(255/np.sqrt(mse)) # 255 is the maximum grey value
return mse, snr, psnr
def convert(img, target_type_min, target_type_max, target_type):
imin = img.min();
imax = img.max();
a = (target_type_max - target_type_min) / (imax - imin);
b = target_type_max - a * imax;
new_img = (a * img + b).astype(target_type);
return new_img
def downsample(sig, factor, Anti_Aliasing=False, img=False):
if Anti_Aliasing:
if img:
lpf = np.ones((3,3))/9
sig = cv.filter2D(sig, -1, lpf)
else:
b, a = butter(15, 1/factor, btype='low', analog=False)
sig = lfilter(b, a, sig)
if img:
sig = sig[::factor, ::factor]
else:
sig = sig[::factor]
return sig
def uint8c(x):
"""custom uint8 - (if val < 0 then val = 0)"""
x2 = np.copy(x)
x2 = x2.clip(0, 255)
return np.uint8(x2)
def sample_test(x,factor, printout=True):
print(f"")
x_downsampled = downsample(x, factor)
x_downsampled_antialiasing = downsample(x, factor, True)
x_fft = np.abs(fft(x))
x_downsampled_fft = np.abs(fft(x_downsampled))
x_downsampled_antialiasing_fft = np.abs(fft(x_downsampled_antialiasing))
f1 = np.arange(-len(x_fft)/2, len(x_fft)/2)
f2 = np.arange(-len(x_downsampled_fft)/2, len(x_downsampled_fft)/2)
f3 = np.arange(-len(x_downsampled_antialiasing_fft)/2, len(x_downsampled_antialiasing_fft)/2)
if printout:
plt.figure(figsize=(15,10))
subplot(321); stem(x); title('original')
subplot(322); stem(f1, fftshift(x_fft)); title('original fft')
subplot(323); stem(x_downsampled); title(f'decimation factor = {factor}')
subplot(324); stem(f2, fftshift(x_downsampled_fft)); title('fft decimation')
subplot(325);stem(x_downsampled_antialiasing);title('decimation antialiasing')
subplot(326); stem(f3, fftshift(x_downsampled_antialiasing_fft)); title('decimation antialiasing fft')
plt.tight_layout()
return x_downsampled, x_downsampled_antialiasing
time = 1 # sec
fs = 200 # Hz
f = 10 # Hz
A = 7
Ts = 1/fs
t = np.arange(0, time-Ts, Ts)
x = A*cos(2*pi*10*t) + 2*A*cos(2*pi*17*t)+3*A*cos(2*pi*5*t);
sample_test(x, 2)
sample_test(x, 4)
sample_test(x, 8)
sample_test(x, 17)
(array([ 42. , -27.76897101, 20.96755674, -17.0690992 ,
-9.52913744, 17.03932492, -21.05361932, 37.90040574,
-27.84310872, 19.41239557, -20.31479123, -2.58612709]),
array([ 5.52411930e-15, 2.57472125e-04, 3.39310392e-01, 4.89487863e+00,
-3.31853335e+00, -6.24926302e+00, 1.44132060e+01, -1.89891894e+01,
1.84453961e+01, -1.38034405e+01, 5.80662156e+00, 3.68393273e+00]))
img = cv.imread(r"C:\Users\rom21\OneDrive\Desktop\imageLAB4\cameraman.tif",0)
x=img[70,:]
sample_test(x, 2)
sample_test(x, 4)
sample_test(x, 8)
sample_test(x, 17)
(array([162, 168, 171, 177, 14, 16, 79, 144, 63, 134, 168, 169, 163,
157, 151, 144], dtype=uint8),
array([ 2.13073173e-14, 1.12114420e-03, 2.08289626e+00, 5.86659563e+01,
1.90954452e+02, 1.66223053e+02, 1.74678454e+02, 9.83228354e+01,
-1.99445170e+01, 6.29735645e+01, 1.10052834e+02, 1.13039900e+02,
1.05583938e+02, 1.91917891e+02, 1.68573857e+02, 1.66273204e+02]))
My signal can see aliasing
x = A*cos(2*pi*10*t) + 2*A*cos(2*pi*20*t) + 2*A*cos(2*pi*30*t)
sample_test(x, 2)
sample_test(x, 4)
sample_test(x, 8)
sample_test(x, 17)
(array([ 35. , -13.52653238, -2.16311896, -3.56015122,
2.98935688, -14. , 5.66311896, 26.21262707,
-20.48935688, 4.87405654, -7. , 4.87405654]),
array([ 4.60343275e-15, 1.62427404e-04, 1.37472032e-01, 1.19119573e+00,
-3.18374634e-01, -2.94432192e-02, 1.19656826e-01, -8.96916238e-02,
7.07075329e-02, -6.54555921e-02, 5.92103586e-02, -3.74673268e-02]))
def upsample1d(sig, factor, method = 'ZeroOrderHold'):
if method == 'ZeroOrderHold':
sig_upsampled = np.repeat(sig, factor)
elif (method == 'FirstOrderHold'):
sig_upsampled = np.zeros(len(sig) * factor)
sig_upsampled[::factor] = sig
filter_linear = np.zeros(factor*2 -1)
for i in range(int((len(filter_linear) + 1 )/2)):
filter_linear[i] = filter_linear[int(len(filter_linear)-1-i)] = (i+1)/factor
#= filter_linear[int(len(filter_linear)-1 -i)]
sig_upsampled = np.convolve(sig_upsampled, filter_linear, 'same')
elif (method == 'sinc'):
sig_upsampled = np.zeros(len(sig) * factor)
sig_upsampled[::factor] = sig
l = 50
x = np.linspace(-l, l, factor*l-1)
sinc_filter =np.sinc((np.pi/(factor+2.1))*x)
sig_upsampled = np.convolve(sig_upsampled, sinc_filter, 'same')
elif (method == 'cubic'):
sig_upsampled = np.zeros(len(sig) * factor)
for i in range(0,len(sig)-3):
f = CubicSpline([0, 1, 2], sig[i:i+3], bc_type='natural')
x_new = np.linspace(0, 2, factor+1)
sig_upsampled[i*factor:i*factor+factor+1] = f(x_new)
return sig_upsampled
x = A*cos(2*pi*10*t) + 2*A*cos(2*pi*17*t)+3*A*cos(2*pi*5*t);
factor = 4
x_downsampled, x_downsampled_antialiasing = sample_test(x, factor, False)
#s3 = upsample1d(x_downsampled, factor, method = 'ZeroOrderHold')
#s4 = upsample1d(x_downsampled_antialiasing, factor, method = 'FirstOrderHold')
#s5 = upsample1d(x_downsampled_antialiasing, factor, method = 'Ideal')
methods= ['ZeroOrderHold', 'FirstOrderHold', 'sinc', 'cubic']
for method in methods:
upsampled = upsample1d(x_downsampled, factor, method = method)
upsampled = upsampled[:len(upsampled)-1]
MSE = np.sum(np.square(x - upsampled) / len(x))
SNR = 10*np.log10(np.std(x) / np.std(x - upsampled))
fig = plt.figure()
axes = fig.subplots(nrows=2, ncols=1)
axes[0].plot(range(0, len(x)), x, range(0, len(x)), upsampled)
axes[0].set(title = 'original signal (blue), upsampled ' + method + ' (orange)\nMSE = ' + "{:.3f}".format(MSE) + ', SNR = ' + "{:.2f}".format(SNR) + 'dB')
axes[1].plot(x_downsampled)
axes[1].set(title = 'downsampled signal')
fig.tight_layout()
def quantize(signal, k):
temp = signal
if type(signal[0]) == np.float64:
temp = uint8c(signal*255)
delta = 256/k
signal_q = uint8c(np.int64(np.double(temp) / delta) * delta + delta/2)
if type(signal[0]) == np.float64:
signal_q = np.double(signal_q / 255)
return signal_q
def quatize_analaze(sig, num_bits):
sig_float_q = quantize(sig, num_bits)
err = sig - sig_float_q
avg_err = np.sqrt(np.average(np.square(err)))
snr_db = 20*np.log10(np.std(sig) / np.std(err))
fig, axi = plt.subplots(2,1)
axi[0].hist(sig.flatten(), bins = 255)
axi[0].hist(sig_float_q.flatten(), bins = 255)
axi[0].set(title = 'Histogram')
axi[0].legend(['Float64','Float64 Quantized ' + str(num_bits) + '-bit'])
axi[1].plot(np.arange(0, len(sig)), sig,np.arange(0,len(sig_float_q)), sig_float_q)
axi[1].set(title = 'Signals\nSNR = ' + "{:.1f}".format(snr_db) + 'dB, Avg Error = ' + "{:.4f}".format(avg_err))
axi[1].legend(['Float64','Float64 Quantized ' + str(num_bits) + '-bit'])
fig.tight_layout()
plt.show()
return sig_float_q
fs = 200
L = 1
f = 2
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
y = uint8c((1 + y) * 128)
Min = np.min(y)
Max = np.max(y)
y_double = ((np.double(y) - Min) / (Max - Min))
img=cv.imread(r"C:\Users\rom21\OneDrive\Desktop\imageLAB4\cameraman.tif", 0)
factor=17
x=img[70,:]
y_4 = y/4 +96
num_bits = 4
quatize_analaze(x, num_bits)
quatize_analaze(y_double, num_bits)
quatize_analaze(y, num_bits)
quatize_analaze(uint8c(y_4), num_bits)
target_noise_db = 10
target_noise_watts = 10 ** (target_noise_db / 10)
mean_noise = 0
noise = np.random.normal(mean_noise, np.sqrt(target_noise_watts), len(y))
noise_sig = y+noise
noise_sig =((np.double(noise_sig) - Min) / (Max - Min))
sig_float_q = quatize_analaze(noise_sig, 2)
a. Downsample factors = [2, 4, 10]:
img = cv.imread(r"C:\Users\rom21\OneDrive\Desktop\imageLAB4\chrono.tif", 0)
chrono = img[1601:2500, 1201:2300]
height, width = chrono.shape
factors = [2, 4, 10]
down_chrono_list = []
#downsampled
plt.figure(figsize=(15,15))
plt.title(f"orig")
plt.imshow(chrono, cmap = 'gray', vmin = 0, vmax = 255, aspect='equal')
for i, factor in enumerate(factors):
downsampled_img = downsample(chrono, factor, img=True)
plt.figure(figsize=(15,15))
plt.title(f"factor:{factor}")
plt.imshow(downsampled_img, cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
For a 1: 2 aspect ratio it will be difficult to see the differences but it will be enough to take a look at the size of the matrix that represents the image to understand that we have shrunk it while losing information along the way of course. For decimation in a ratio of 1: 4 and of course for decimation in a ratio of 1:10 we can see the effects of matrix contraction.
def Zero_Order_Hold2d(img, factor):
new_img = np.zeros(tuple(np.array(img.shape)*factor))
for c in range(img.shape[0]):
new_img[c*factor:c*factor+factor,:] = np.repeat(img[c,:], factor)
return new_img
def bilinear(image, factor):
img_height, img_width = image.shape
height, width = tuple(np.array(image.shape)*factor)
image = image.ravel()
x_ratio = float(img_width - 1) / (width - 1)
y_ratio = float(img_height - 1) / (height - 1)
y, x = np.divmod(np.arange(height * width), width)
x_l = np.floor(x_ratio * x).astype('int32')
y_l = np.floor(y_ratio * y).astype('int32')
x_h = np.ceil(x_ratio * x).astype('int32')
y_h = np.ceil(y_ratio * y).astype('int32')
x_weight = (x_ratio * x) - x_l
y_weight = (y_ratio * y) - y_l
a = image[y_l * img_width + x_l]
b = image[y_l * img_width + x_h]
c = image[y_h * img_width + x_l]
d = image[y_h * img_width + x_h]
resized = a * (1 - x_weight) * (1 - y_weight) + \
b * x_weight * (1 - y_weight) + \
c * y_weight * (1 - x_weight) + \
d * x_weight * y_weight
return resized.reshape(height, width)
#%% 2.1.3 Interpolation Zero_Order_Hold 2d
for i, factor in enumerate(factors):
downsampled_img = downsample(chrono, factor, img=True)
resized = Zero_Order_Hold2d(downsampled_img, factor)
#resized = cv.resize(downsampled_img, img.shape, interpolation=cv.INTER_LINEAR)
plt.figure(figsize=(15,15))
plt.title(f"factor {factor}, Zero_Order_Hold2d Interpolation")
plt.imshow(resized, cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
#%% 2.1.3 Interpolation Bilinear 2d
dpi = [1250, 300, 150, 72]
for i, factor in enumerate(factors):
downsampled_img = downsample(chrono, factor, img=True)
resized = bilinear(downsampled_img, factor)
plt.figure(figsize=(15,15))
plt.title(f" Bilinear Interpolation {factor}")
plt.imshow(resized, cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
image_path = r"C:\Users\rom21\OneDrive\Desktop\imageLAB4\skull.JPG"
img = cv.imread(image_path)
img_quarnt=[]
quant_levels = [2, 4, 8, 16, 32, 64, 128]
for q_level in quant_levels:
img_quarnt.append(quantize(img, q_level))
fig,ax = plt.subplots(2,4,figsize=(15,15))
ax[0,0].set_title("orig")
ax[0,1].set_title("quant level = 128")
ax[0,2].set_title("quant level = 64")
ax[0,3].set_title("quant level = 32")
ax[1,0].set_title("quant level = 16")
ax[1,1].set_title("quant level = 8")
ax[1,2].set_title("quant level = 4")
ax[1,3].set_title("quant level = 2")
ax[0,0].imshow(img, cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[0,1].imshow(img_quarnt[6], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[0,2].imshow(img_quarnt[5], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[0,3].imshow(img_quarnt[4], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[1,0].imshow(img_quarnt[3], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[1,1].imshow(img_quarnt[2], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[1,2].imshow(img_quarnt[1], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[1,3].imshow(img_quarnt[0], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
plt.tight_layout()
plt.show()
#%%
fig,ax = plt.subplots(2,4,figsize=(15,10))
ax[0,0].set_title("orig")
ax[0,1].set_title("quant level = 128")
ax[0,2].set_title("quant level = 64")
ax[0,3].set_title("quant level = 32")
ax[1,0].set_title("quant level = 16")
ax[1,1].set_title("quant level = 8")
ax[1,2].set_title("quant level = 4")
ax[1,3].set_title("quant level = 2")
ax[0,0].hist(img.flatten(), 256, [0, 256])
ax[0,1].hist(img_quarnt[6].flatten(), 256, [0, 256])
ax[0,2].hist(img_quarnt[5].flatten(), 256, [0, 256])
ax[0,3].hist(img_quarnt[4].flatten(), 256, [0, 256])
ax[1,0].hist(img_quarnt[3].flatten(), 256, [0, 256])
ax[1,1].hist(img_quarnt[2].flatten(), 256, [0, 256])
ax[1,2].hist(img_quarnt[1].flatten(), 256, [0, 256])
ax[1,3].hist(img_quarnt[0].flatten(), 256, [0, 256])
fig.tight_layout()
plt.show()